Magnetostatics solution (nonlinear case)ΒΆ

[1]:
%%capture
%run current_density.ipynb
%run TEAM13_nonlinearity.ipynb
[2]:
##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)
[3]:
R
[3]:
<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
        with 23597 stored elements in Compressed Sparse Column format>
[4]:
##############################################################################
# Assembly
##############################################################################

linear = '*coil,default'
nonlinear = 'r_steel,l_steel,mid_steel'
maxIter = 100

order = 1

fem_linear = pde.int.evaluate3(MESH, order = order, regions = linear).diagonal()
fem_nonlinear = pde.int.evaluate3(MESH, order = order, regions = nonlinear).diagonal()

D = pde.int.assemble3(MESH, order = order)

phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)

aJ = jx_L2 @ D @ phix_Hcurl.T +\
     jy_L2 @ D @ phiy_Hcurl.T +\
     jz_L2 @ D @ phiz_Hcurl.T

aJ = 1e7*aJ

def gss(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A

    Kxx = curlphix_Hcurl @ D @ sp.diags(fxx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
    Kyy = curlphiy_Hcurl @ D @ sp.diags(fyy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
    Kzz = curlphiz_Hcurl @ D @ sp.diags(fzz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T

    Kxy = curlphiy_Hcurl @ D @ sp.diags(fxy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
    Kxz = curlphiz_Hcurl @ D @ sp.diags(fxz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T

    Kyx = curlphix_Hcurl @ D @ sp.diags(fyx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
    Kyz = curlphiz_Hcurl @ D @ sp.diags(fyz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T

    Kzx = curlphix_Hcurl @ D @ sp.diags(fzx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T
    Kzy = curlphiy_Hcurl @ D @ sp.diags(fzy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T

    return Kxx + Kyy + Kzz + Kxy + Kxz + Kyx + Kyz + Kzx + Kzy

def gs(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
    return curlphix_Hcurl @ D @ (fx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
           curlphiy_Hcurl @ D @ (fy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
           curlphiz_Hcurl @ D @ (fz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) - aJ

def J(A):
    curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
    return np.ones(D.size)@ D @(f_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + f_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) -aJ@A




import time
A = np.zeros(curlphix_Hcurl.shape[0])
mu = 0.0001
eps_newton = 1e-5
factor_residual = 1/2

tm = time.monotonic()
for i in range(maxIter):

    gssu = R.T @ gss(A) @ R
    gsu = R.T @ gs(A)

    tm = time.monotonic()
    # wS = chol(gssu).solve_A(-gsu)
    wS = sp.linalg.spsolve(gssu, -gsu)
    print('Solving took ', time.monotonic()-tm)

    w = R@wS

    alpha = 1

    # ResidualLineSearch
    # for k in range(1000):
    #     if np.linalg.norm(gs(A+alpha*w),np.inf) <= np.linalg.norm(gs(A),np.inf): break
    #     else: alpha = alpha*factor_residual

    # AmijoBacktracking
    float_eps = 1e-12; #float_eps = np.finfo(float).eps
    for kk in range(1000):
        if J(A+alpha*w)-J(A) <= alpha*mu*(gsu@wS) + np.abs(J(A))*float_eps: break
        else: alpha = alpha*factor_residual


    A_old_i = A
    A = A + alpha*w

    print ("NEWTON: Iteration: %2d " %(i+1)+"||obj: %2e" %J(A)+"|| ||grad||: %2e" %np.linalg.norm(R.T @ gs(A),np.inf)+"||alpha: %2e" % (alpha))


    # if ( np.linalg.norm(R.T @ gs(A),np.inf) < eps_newton):
    #     break
    if (np.abs(J(A)-J(A_old_i)) < 1e-5):
        break

elapsed = time.monotonic()-tm
print('Solving took ', elapsed, 'seconds')
Solving took  2.796000000089407
NEWTON: Iteration:  1 ||obj: -3.220290e+10|| ||grad||: 2.479687e+07||alpha: 6.250000e-02
Solving took  2.5620000001508743
NEWTON: Iteration:  2 ||obj: -6.401723e+10|| ||grad||: 5.435596e+07||alpha: 1.000000e+00
Solving took  2.672000000020489
NEWTON: Iteration:  3 ||obj: -6.639095e+10|| ||grad||: 5.773030e+07||alpha: 5.000000e-01
Solving took  2.6410000000614673
NEWTON: Iteration:  4 ||obj: -1.508296e+11|| ||grad||: 6.567249e+07||alpha: 1.000000e+00
Solving took  2.5780000002123415
NEWTON: Iteration:  5 ||obj: -1.530142e+11|| ||grad||: 4.925437e+07||alpha: 2.500000e-01
Solving took  2.547000000020489
NEWTON: Iteration:  6 ||obj: -1.858670e+11|| ||grad||: 3.979715e+07||alpha: 1.000000e+00
Solving took  2.6089999999385327
NEWTON: Iteration:  7 ||obj: -1.878539e+11|| ||grad||: 3.481956e+07||alpha: 1.250000e-01
Solving took  2.515000000130385
NEWTON: Iteration:  8 ||obj: -1.909770e+11|| ||grad||: 2.171232e+07||alpha: 2.500000e-01
Solving took  2.594000000040978
NEWTON: Iteration:  9 ||obj: -1.931830e+11|| ||grad||: 1.881516e+07||alpha: 5.000000e-01
Solving took  2.6089999999385327
NEWTON: Iteration: 10 ||obj: -1.933408e+11|| ||grad||: 2.100793e+07||alpha: 2.500000e-01
Solving took  2.5310000001918525
NEWTON: Iteration: 11 ||obj: -1.936156e+11|| ||grad||: 1.285842e+07||alpha: 1.000000e+00
Solving took  2.5
NEWTON: Iteration: 12 ||obj: -1.939955e+11|| ||grad||: 6.739643e+06||alpha: 2.500000e-01
Solving took  2.547000000020489
NEWTON: Iteration: 13 ||obj: -1.940739e+11|| ||grad||: 3.975133e+06||alpha: 5.000000e-01
Solving took  2.530999999959022
NEWTON: Iteration: 14 ||obj: -1.941111e+11|| ||grad||: 1.127736e+06||alpha: 1.000000e+00
Solving took  2.5
NEWTON: Iteration: 15 ||obj: -1.941138e+11|| ||grad||: 1.987664e+05||alpha: 1.000000e+00
Solving took  2.5
NEWTON: Iteration: 16 ||obj: -1.941138e+11|| ||grad||: 9.657019e+03||alpha: 1.000000e+00
Solving took  2.530999999959022
NEWTON: Iteration: 17 ||obj: -1.941138e+11|| ||grad||: 7.475946e+00||alpha: 1.000000e+00
Solving took  2.5469999997876585
NEWTON: Iteration: 18 ||obj: -1.941138e+11|| ||grad||: 1.351817e-02||alpha: 1.000000e+00
Solving took  2.547000000020489
NEWTON: Iteration: 19 ||obj: -1.941138e+11|| ||grad||: 2.625993e-05||alpha: 1.000000e+00
Solving took  2.5939999998081475 seconds
[5]:
curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
Bx = curlphix_Hcurl_P0.T @ A
By = curlphiy_Hcurl_P0.T @ A
Bz = curlphiz_Hcurl_P0.T @ A
[6]:
##############################################################################
# Storing to vtk
##############################################################################

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')

[7]:
import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,4])

p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)

threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=18)
p.add_mesh(arrows, color="orange")


mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")

p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")
[ ]: